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ABSTRACT 


One case of analytic data and one case of real data were numerically 
integrated using a 5-level baroclinic primitive equation global model of 
the general circulation. The feasibility of using a mixed second order 
and fourth order space differencing scheme to improve the phase speeds 
of meteorological waves was examined. The results indicate that mixed 
scheme tends to give a better representation of the phase speeds than the 


second order scheme does. 
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I. INTRODUCTION 


A baroclinic primitive equation model using a global, staggered, 
spherical, sigma coordinate system which was originally written by 
Dr. F. J. Winninghoff and which is now being developed by the United 
States Navy Fleet Numerical Weather Central (Elias, 1973) was used in 
this study. The purpose of this investigation was to modify the hori- 
zontal space difference equations to a mixed second and fourth order 
scheme with the aim of improving the phase speeds of the meteorologi- 
cal waves, as proposed by Williams (1972). Real time data (Elias, 
1973) and analytic data using an analytic spherical harmonic stream 
function (Neatam, 1946) were used as initial conditions. The use of 
the analytic case allowed the control of temperature and moisture dis- 


tributions, predominant wave number, phase speed and wave amplitude. 


del. 





II. BAROCLINIC PRIMITIVE EQUATION MODEL 


The governing differential equations, written in vector form, are 
similar to sets used by Smagorinsky et al (1965), Arakawa, et al (1969) 
and Kesel and Winninghoff (1972). The integrations were carried out on 
a global, spherical, staggered, sigma coordinate system using the con- 
servative-type difference equations based on the procedure of Arakawa 
(1966). The complete set of finite difference equations are given in 
Kesel and Winninghoff (1972) for the non-staggered, polar stereographic 
grid. Appendix A contains the finite difference equations for the 
spherical, sigma, staggered system of equations. The heating, moisture, 
friction and diffusion finite difference terms are given in Kesel and 


Winninghoff (1972). 


A. PRIMITIVE EQUATIONS 
The following are the equations for the model in spherical, sigma 
(o) coordinates. 


The dimensionless vertical coordinate is: 
Gi= P/l . 


The vertical velocity measure is defined as: 


The latitudinal component of the equation of motion is: 


out a 1 jo(uuT) a ¢ Grins 
ot a cos 9 OA ele) 
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The longitudinal-component of the equation of motion is: 


OvT _ 1 potuvT) 8 Gmieose 
t a cos 8 OA ele) 


. Td (wv) fe Tuutand _ 


7G Tuf 


1 ,.7 OO , RT OT 
Bellis act a ee ee 


The First Law is: 


ont 1 0 (TuT) d(TvTcos8) 











ont _ _ ee 

ot a cos 6 OA ag 00 
A To (wT) 

00 

RT oT 1 u oT om cos 8 
* e0 ee eeeecos cd) ON 30 })] 
+ TH + De. ; 

The Moisture Equation is: 
oqT _ 1 pole) i d(Tvqcos®) | 
ot a cos 6 OA 00 
T 9(wq) 

+ 5 On + ae ° 
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(1) 


(2) 


(3) 


(4) 





The continuity equation is: 


2 an 1 d(uT) , d(vtcosé@) T OW 

hae acc ie. m weoo 9 ~=(Sc - (5) 
The Hydrostatic Equation is: 

od _ RT 

i tamer (6) 


Fj 


In the above equations, u is the zonal wind component; v, the 
meridional wind component; 1, the terrain pressure, T, the temperature; 
f, the Coriolis parameter; ¢%, the geopotential; q, the specific humidity; 
H, the diabatic heating; F, the frictional stress; D, the lateral diffu- 
Sion of momentum; Drs the lateral diffusion of.heat; nee the lateral 
diffusion of moisture source/sink term; P, the pressure; a, the radius 
of the earth; A, longitude and 0, the latitude. 

Kesel and Winninghoff (1972) modified the pressure gradient terms 
as suggested by Kurihara (1968) to reduce the inconsistent truncation 
error which arises as a result of the sensitivity of the model to varia- 
tions in the height and character of the terrain. The value of the geo- 
potential is interpolated at neighboring gridpoint locations to the same 
pressure surface as at the computation point to enable one to compute 
to geopotential gradients on pressure surfaces which are thus synthe- 
Sized locally around each computational point. The pressure gradient 


term in sigma coordinates is: 


4 RT 
me seo vi (7) 
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Solving for the gradient of 7 in component form one gets, 


on 71 do dd 

5, ~ pr bx - 37 ] (8) 
Pp re) 

on om 8 _ 36 

20 RT '30 70! (9) 
P oO 


It is shown in Appendix A that the right hand side of equations (8) and 
(9) are the sum of the thicknesses between the sigma surface and the 
pressure surface to the right and left of the computational points. 
Since the thickness depends only upon the mean temperature of the layer 
in question, the thickness equation is used to solve equations (8) and 
(9). The mean temperature is obtained by interpolation as indicated in 


Appendix A. 


B. GRID 

The longitudinal and latitudinal grid increments were ten degrees 
for the mass and velocity fields. The geopotential (6), temperature 
(T), specific humidity (q), and the vertical-velocity (w = -G) were 
carried at the poles. When 9%, T, q, w and 7 were needed at a velocity 
point they were taken as the average of the four surrounding mass points. 


Figure 1 represents the staggered grid and location of the variables. 


C. VERTICAL LAYERING 

Figure 2 is the vertical layering as described by Elias (1973). 
The basic variables, as indicated, are carried at the center of each 
layer except that the specific humidity (q) is carried at the lower three 


layers only. The vertical velocity (w) is carried at the interface 
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between layers. The vertical coordinate, sigma (0) was defined by 


Phillips (1957). 


D. TIME DIFFERENCING 
The two. finite difference techniques that were used were the cen- 

tered difference scheme and the Matsuno (Euler-backward) difference 
scheme. The finite difference equation for the centered technique is 
t+1 t-l OF 


F = F + 2At —— 
ot 


and that for the Matsuno scheme is 


t 
reo po 4 vt OF 
ot 
ie 

pitl _ poy ae OF 
ot 


Since the centered scheme is not feasible for the first time step and 
also produces a computational mode (Haltiner, 1971), the Matsuno scheme 
was used for the initial time step of each 3~-hour integration to reduce 
solution separation. The Matsuno scheme also selectively dampens high 
frequency waves (Haltiner, 19719 The Arakawa technique (Langlois and 
Kwok, 1969) of averaging quantities involved in the longitudinal deri- 
vitives is used. This allowed a time step of ten minutes to be used. 
Without this technique, the Von Neuman linear computational stability 
criterion (Haltiner, 1971) would require a two and one half minute time 
step because of the longitudinal grid distance at 85° latitude which is 


only 47 kilometers. 
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The lateral diffusion, friction, convective condensation and large 
scale condensation are computed at every time step while heating is com- 


puted every six time steps (Kesel and Winninghoff, 1972). 


E. FOURTH ORDER MIXED SPACE DIFFERENCING 

Williams (1972) found that the stability criterion for fourth order 
unstaggered space differencing required a time step which was 27% smaller 
than that for second order space differencing. He then noted that the 
phase speed of the meteorological wave is independent of the difference 
approximation which is used for the pressure force term, A mixture of 
fourth order differencing in the advection terms and the divergence term 
in the continuity equation and second order differencing in the pressure 
force term was proposed. The stability criterion for this mixed scheme 
was found to require a time step that was only 14% smaller than that of 
the second order scheme. Thus the mixed differencing allows a larger 
time step while giving fourth order accuracy in the phase speed of the 
meteorological wave. Appendix B presents the fourth order differencing 
of the advection terms and the divergence term of the continuity equa- 


tion. A time step of 8 minutes was used for this scheme. 


19 





De -—° + = 












—ee =. _ ~4> a=. 
= =e G ME = 
= - @ =e bbe 








III. INITIAL CONDITIONS 


Two cases, one using real data and the other analytic data, were 
used for initialization of both the second order and mixed differencing 
schemes. The real data was produced external to the main program on a 
Northern Hemispheric FNWC 63 x 63 grid. The analytic data with known 
phase propagation was produced following the work of Neamton (1946), 
Gates (1962), Heburn (1972) and Elias (1973). The real data were fields 
analyzed by FNWC objective schemes. 

The main program was used to interpolate the data so that values as 
5° latitude-longitude intersections were available on the Northern Hemi- 
sphere (Elias, 1973). The data was then reflected into the Southern 
Hemisphere. In both cases, the input data consisted of temperature 
analyses for the Northern Hemisphere at 12 constant pressure levels from 
1000 to 50 mb, height analyses at ten of these levels, moisture analyses 
at four levels from the surface ys 500 mb, sea level pressure and sea 
surface temperature. Monthly mean surface temperature fields were used 
to derive an albedo field (Dickson and Posey, 1967). In both cases the 
terrain height was set to zero. For the barotropic analytic case heat- 
ing, convective adjustment and friction were turned off. 

If one assumes nondivergent horizontal flow, harmonic wave solutions 
of the complete vorticity equation can be obtained for the sphere 
(Neamton, 1946). These solutions give the velocity of propagation of 


the wave. The solution of the ~ field (Elias, 1973) was found to be 


S = A sin (m\-vt)P™ (sin 8) - B sin 6+ C Po (sin 6) (1) 
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where A, B and C are constants to be determined, a is the radius of the 
earth, m is the hemispheric wave number, v/m is the angular phase speed 
of the wave (radians/second), PO denotes the Legendre polynomial and Ba 
represents the associated Legendre function. 


If C = 0 and n = m+ 1 (Elias, 1973), equation (1) reduces to 


YW =A sin (md, - vt) (2N!/25N!)sin 8 (cos 9)° erg sin 6 (2) 


The constant A is arbitrary and proportional to the wave amplitude. It 
was shown by Haurwitz (1946) that the solution obtained for implies 
the existence of a velocity distribution over a sphere such that the 
angular velocity of the westerly current is made up of the sum of three 
terms. The first, B, is constant over the sphere; the second term 
varies along a meridian but is constant around a latitude circle; and 


the third term represents the angular phase speed which is given by 


n(ntl)=-2 _ 22 


vires = & n(ntl) = n(nt1) 


(3) 


where 2 is the earth's rotation. For the analytic experiment A was chosen 
2 -l -/ =] : 
as 1000 m sec } m as 6, n as 7 and B as 8 x 10 sec. The solution of 
equation (3) gives an angular velocity of -9° longitude/day (toward the 
west). Since wave number six was used, integrations were performed for 
only one-sixth of a latitude circle with the remaining five-sixths being 
set equal to the first part. This resulted in a 40% saving of computer 
time for the analytic case only. 


The geopotential fields were derived by solution of the linear 


balance equation and the temperatures, constant at each pressure surface, 


fa 





were consistent with the National Advisory Committee for Aeronautics 
(NACA). Both the geopotential and temperature fields were derived as 
shown by Elias (1973). 

The wind field was produced following the work of Elias (1973) by 


using the finite difference expressions 


_ 1a 

u = a AO (4) 
__ 1 & 

aT (5) 


where a is the earth's radius, A® is the latitudinal distance increment 
and AA is the longitudinal distance increment. Poleward of 25° latitude 
the initial geopotential field, ?, was retained and the -field was 


obtained by the solution of the linear balance equation, 


Vy + Vy > VE/E = V6 /8 (6) 


Equatorward of 25° latitude the ~-field was determined using, 


y=  /f 


where f is a mean Coriolis parameter. A new geopotential field (9%) 


was obtained using 


2 2 
Vo= (£V P+ VW ° VE). 


At 25° latitude, wW and $, were a combination of half the value of the 


poleward case and half of the equatorward case. 
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IV. WAVE ANALYSIS METHOD 


A Fourier series was determined at each five degrees of latitude 
around the latitude circle. With this technique, the phase angles and 
amplitudes of each wave number around a latitude circle can be calcu- 


lated. A Fourier series can be expressed as follows (Heburn, 1972): 


F(x) = A, + > (A, cos mx + BF sin mx) 
m 
= C, +>) (Cc. cos (mx — 6m)) 
m 
on . 
waste eis ay em) 
and 
-1 
je ean (B/A). 
m nm m 


The analytic experiment involved an input stream function of wave 
number six. The values of primary interest were 56 and Ce which are 


the phase angle and amplitude of wave number six. Other angles and 


amplitudes were extracted, especially in the real data case. 
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V. RESULTS 


Two experiments will be presented in this section. The first 
experiment was performed with analytic fields and the second with FNWC 
analyzed fields from 1200Z 10 May 1973. 

As indicated in section III the angular phase speed for wave number 
six was -9° longitude/day (toward the west). A 72-hour forecast of 
surface pressure was made by both the second and mixed order differenc- 
ing. The changes in phase angle vs latitude were plotted for both 
cases and are shown in Figures 3 and 4. Table I summarizes the average 
angular movement per 24-hour the entire 72-hour forecast for both dif- 
ferencing schemes. It should be noted that the mixed differencing 
scheme propagated the wave slower to the west than the second order 
scheme. The barotropic Rossby equation (Holton, 1972) is given as 


B 


ee 
fe 


where £8 is a f is the Coriolis parameter, u is the mean zonal flow, 
k is the zonal wave length, m is the wave number and C is the zonal 
phase speed. Since the quantities 8, k and m are constant the phase 
speed, C, depends only upon the mean zonal flow. Fourth-order differ- 
encing should give a better approximation for the effect of the mean 
zonal wind, u, and therefore the westward movement of the wave would 


be slower than that of the second order differencing case. It would 


appear that the mixed differencing scheme does give a better solution. 
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Phase Angle vs Latitude for the Second Order Scheme. 
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SECOND ORDER MIXED SCHEME 


LATITUDE DEG LAT/DAY DEG LAT/DAY 
70° -12.6° -4.8 
60° -10.3° -7.9 
50° -9,2° -7.4 
40° -8.5° -/7.5 
30° -7.7° -/.2 
20° -7.3° | -7.6 
10° -7.0° -6.5 


0° -4.6° -3.6 


TABLE I. Average Movement of Analytic Wave per 24-Hours for a 
72-Hour Forecast for Both Second Order Differencing and 
the Mixed Scheme. 
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A 48-hour forecast was also made by both the second order and 
mixed differencing schemes using analyzed FNWC data fields. Chart A 
is the initial analyzed surface pressure field, Charts B through E 
are the 12-, 24-, 36- and 48-hour forecasts produced by the second 
order scheme and Charts F through I are the corresponding forecasts 
produced by the mixed differencing scheme. Table II lists five high 
pressure centers and Table III lists seven low pressure centers with 
their initial position and both the 48-hour second order and 48-hour 
mixed scheme forecast positions. As can be seen each high pressure 
center was forecast at least five degrees farther to the east by the 
mixed differencing scheme. Low pressure centers 1 and 3 also were 
forecast five degrees farther to the east by the mixed scheme while 
low centers 2 and 4 were moved very little by both differencing schemes. 
Low number 5 was the only low that was moved slower than the second 
order scheme. Low number 6 was developed into two centers by both 
forecast schemes. Only the forecast position of the major center was 
listed in Table III. The mixed scheme moved this center nine degrees 
farther than the second order scheme. Center number 7 was also 
initially a single center but was developed into three small centers 
with the second order forecast scheme. The mixed scheme maintained 


the more realistic single center system near its initial position. 
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36-Hour Forecast, Second Order Scheme. 
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Chart E. 
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\i 


ta 
Ww 


wat | 0.oO oO 


PROJECTION POLAA STLILOGRAPHIC-TRUL AT a” WORTH 





36 


36-Hour Forecast, Mixed Scheme. 
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HIGH CENTER INITIAL 48-HOUR SECOND 48-HOUR MIXED 


NUMBER POSITION ORDER POSITION SCHEME POSITION 
“4 37N 177W 37N 177W 35N 172W 

2 35N 136W 47N 128W 43N 121W 

3 None 31N 53W 33N 48W 

4 45N  25W 52N 17W 52N 12W 

5 64N 50E 60N 60E 52N 66W 


TABLE II. Initial Positions of High Centers and Forecast Positions 
of the Second Order Scheme and the Mixed Scheme. 


LOW CENTER INITIAL 48-HOUR SECOND 48-HOUR MIXED 
NUMBER POSITION ORDER POSITION SCHEME POSITION 

1 47N 155W 48N 152W 49N 147W 

2 47N 89W 5O0N 45W 48N 40W 

3 47N SOW 50N 87W 51N 88W 

4 67N O2E 68N O3E 68N O4E 

) 61N 89E SON 99E 61N Q97E 

6 36N 152E SON 151E 48N 160E 

7 30N 80E * 33N 81E 


TABLE III. Initial Positions of Low Centers and Forecast Positions 
of the Second Order Scheme and the Mixed Scheme. 


*&Several small lows were forecast by the second order scheme. 
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VI. CONCLUSIONS 


Both the analytic and real data experiments indicate that the 
mixed fourth and second order differencing scheme gives a more 
accurate solution for the meteorological phase speeds. More analytic 
experiments with a more realistic phase speed will have to be per- 
formed to produce the final conclusive proof of the mixed scheme. 

A series of real data experiments should also be considered. In 
addition the fourth order differencing must eventually be extended 


to the poles. 
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APPENDIX A 


DIFFERENCE EQUATIONS 


The longitudinal momentum equation terms in difference form are: 


n+l n-1 
Tu,. — u., 
gut = 14 ia (1) 
ot 2At ¢ 
y. i polit + dtuvT) cost, _ 
a cos @ or els) 
: an at 1,pe" T41/2,4 7 ~ (yu, ij) “a a pun "4H 1254 
a cos 8, 2Ar 
‘i (u, nat ay uJ yavT i eee O 547 - (us 57 ae e _1 C08 ae = @) 
where 
265) = mau ist Tis, sien, 9? Ti, 541%a, 47 "a,5-1%s,5-D ° 
(251) 
ut =1Gg u,.t+T7 u +7 u 
1 2 uadigi | 1—1, jo ee le lt 
Tl fem nee pa a) 
ut =i u..+ 7 
10a fa al B 4° ij ij i,jt+2™ i,jt2° Tea jai" eet 


i,j-1"4,j-V ona 
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1 
od Lae 
ono 
=I 
¢ 
~ 
=I 


Ee? fej? 1, j-1 4-1, j-1 * 


Ty page eer (2.4) 


The longitudinal grid increment (AA) and the latitudinal grid increment 


(A8) are both ten degrees. 


, owu . = age Mi etd 2 Cee Ce ey (3) 
96 ij 2Ao 


where the vertical velocity, -6 » is defined as, 
e eal N 
W = 
k+1 Ge ij 


m 1 pith ik ig gk ijk A THIET Ba 





a cos 6, Ar Aé 
a (4) 

and the local change as 

QMN 25> | OM ag Mayk n VMAj+TKCOSY j 4a Vj 14099951, , 

Ge ij >> Ada cos 8, A@ a cos o. (5) 
The um terms of Equations (4) and (5) are, 
UT. .4. = UL,4.(7,,+ 7 =P oF ae ill ee Coa) 

i+1j i+ ij it+lj ij+l ij-1’ 4 ; 
The quantities Gee uti and uT. are also determined in the 


tli’ ij+l ij-l 


same manner, 
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Tvf = 2 Ls ios (6) 


where 
i. a +1 ) (6.1) 
ij 4“"ag’ "i-1g7 “i-1,j+) “i-1,35-1 
Tuv tanO _ u, (1, .) tand. 
eee (7) 
a a 
] RT ov , 7 OO, _ 
~ a cos 0 pe Bre - 
d..- 6, .) 
1 - ,om* = ( ij i-1, j 
~ a cos e, oo ay ua (8) 
where 
T.. = Lig +T..+T 2 ) (8.1) 
a) OG ad-1j. ig) i-1,4+1 ~“i-1,j-1 
TT = 45 +1.. +7. 4 .,,¢+T7 ) (B52) 
43 ~ Gig” "ig * "i-2,541 7 “i-1,J5-1 : 


From Equation (9) of section II the change of terrain height is found 


as follows 


(6-0) (®&-6) 
om x T R i R L 
Gy = ee ag (9) 
(® - O )+ (6 - @ ) 
See eo an 
RT Ax 
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where %, is to the right of the computational point and e is to the 
left of the computational point. Sigma and p indicate the surface in 
question. The thickness between the pressure surface and the sigma 


surface to the right and left of the computational point is defined as, 


Ao, = (, ~ %, ) and (11) 
Op P Oo 

Ao = (> - 9% ) (12) 
Lop Ly L, 


therefore one can write, 


= 1 

Ad = -RT &n — (13) 
Sap Op P» 
a 5 

Ad, - RT, Ln — (14) 
Op P Po 


The mean temperature is found as follows, 


— dL aekn 1 
~ “comp oo 2 (15) 

and in difference form 

(Tp - Tt) (LN T.. - LN T..) 
T =-7T + Pere te 1) (16) 
oe ij Ao 2 

(Ty - TT) (LN T.. - LN T,..) 
i 7 ‘ o+1 Oo ij ij (17) 
L ij Ao 2 

Op 
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7 =+r. + 


Meg ety 8 1-1,441 * "4-1, 5-17 (17.1) 


IN T,. = =(2n ue + Wn 7. + on Tl , = ¢ er Cl yZ) 


ij i-1lj i-1, j+l1 Agee 
Figure 5 illustrates the differencing in both the horizontal and verti- 
cal, Substituting Equations (16) and (17) into Equations (13) and (14) 


one obtains, 








eee Re (2n 7, .-£n 13) 
—- pir ey ee ee) a i ee 
Ag, RIT, + AG 5 ]-(2n M4 Qn ™5 (18) 
Op 
_ Game i. Conese Tee _ 
Ag, = ey a 5 ] RUT Fae (19) 
Op 
where 
Pp (oe TD a 
dn, 1 = ’n —=——" = Ln Te. = Ln Ts. (20) 
Potl +1 : : 
and 
Pp 
Qn, — = int, - int, _,, (21) 
Potl ! 
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Toome 


—_— =aeSP SE See 





Interpolation for Vn. 


Figure 5, 
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The latitudinal-momentum equation terms in difference form are, 


N+1 N-1l 
avn ag GG a 
ta 2At 
= 1 po drvw) a o(ivvcos®) | = 
a cos @ OX ele) 7 
eal rag its MT i41/25 gi? ws 
a cos = 2AA 
(v;.tv,..5)vT,.,,c088. 15 = (v,.-V.._5)vM,,_,cos8,_) 
i a eS (2) 


The ull terms are determined in the same manner as the terms for 


Equation (2) were found. 


Ta(wv)  _ re wo ajo41 Vij? _ on-1 Vaio iJ Yijo-), (3) 
d0 eee 2Ao 2A0 


The vertical velocity, w, and the local change, au Fe » are obtained 


as in Equations (4) and (5) of the longitudinal equation of momentum. 


tf =f. T.. v.. (4) 


Ts = Equation (6.1) or the u-equation of momentum + 


= u,.(T..u,,)tand, 
Tuu tanO _ i ( ij i ) 


a a 


(5) 
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1 ,w od , RT OT, _ 
se! 76 + 76) = 
(> re 

1 = ,oT,* | = i-1, j+1 i-1,j-1 

"a cos 8 (RTS 5G * 45 ne ke? 
(@ -@ )+(O -& ) 

(aT * - ee (7) 
307, RT Ae 


Equation (7) is solved by the same method as was shown for the 


lontitudinal momentum equation. 


Thermodynamic Equation 





N+1 N-1 
ant) — “Tajo 7 “Tigo ss 
ot " 2At 
a (out) i d(tvTcos8 } = 
a cos 9 oA a0 a 
1 Jit j Minny "iy? - Mga tae 
a cos a, Ar 


T cos, 


Meprigmeerere come), 1 maged 14-4 aq-1°°° "5 nen 
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where T is the temperature at the velocity point based on the average 
of the four surrounding mass points and 7 is also calculated in the 


Same manner. 
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7 O(wI) _ 7 Mijo*ijo41 " Ti 40° 7 i jo-1 “140 i Yijo-)? , (3) 
30 ij 2Ao0 2AG 
. om, N 

The vertical velocity, w, and the local change, GP y » are obtained 
as in the u-equation of motion. 
RT OTT 1 oT OT = 
Cp [wo GE + a cos alt ay t 30 cos0)) | 2 
RT (we pa’ -1) eG (27) 
Cho 2 ij Oo ‘dt’ ij 
_%o Miva Tey 7 Yai ay? 
a cos 6. 2AA 

J 

(V5 oT 44 6089.4, — V5 474410089. _1) 

20 (4) 


where the vertical velocity and the local change are again computed 


as before. 


The Moisture Equation 


ql... - qu.. 
dam = ij (1) 
dt 2At 
pee epoeld .. otvigeos®), - 
a cos 8 oA a0 7 
1 Aisa Missy 2 Shiga lage P 
a cos oe Ar 
Pe etal itl 1a joa Vj 4j-1°°" 5-1, 
Ao 


(2) 
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where q is the specific humidity at the velocity point based on the 
average of the four surrounding mass points and 7 is also calculated 


in the same manner. 


T O(wq) _ 7 ajo Fijot1 a qi jo? 
30 ij 2A0 
_ Mijo-iijo * Yijo-1?, (3) 
2Aoq 
The Hydrostatic Equation 
dd = ~- RT din a 








O 0. 

k kt] 
er a ial = a 
k+1 k O. no, + Ano, 44 — 22nT 


At the lowest level % is given as 


os = - RTs sk Qn (0.9) 
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APPENDIX B 


In accordance with the work of Williams (1972) the fourth order 
terms are a sum of a weighted average of the second order and fourth 
order differencing. The weighting factor the second order part is 
four-thirds and a minus one-third for the fourth order portion. The 


advection term of the longitudinal momentum equation is: 


== 


: 1 4 Mast O54 42UT 41/2: - (u, + U5 _14UT; 105 
a cos oe 3 2AA 
: (u, + U..45)cos8. VT way -(u, .- us,_5)cos0, | wTijel, 
2AA 
: esl Us. 5.)uT. 1. Cu, + U,_5,)UT. 1, 
3 4Ar 
, (u, + u,.4,)c0s9. VT, a5 —(u, + u,._,)cos®,_. Ti i-2) a 


4A 


The vm terms are determined as shown in Appendix A. The advection term 
for the latitudinal momentum equation is similar to term (1). The 


advection term of the thermodynamic equation is: 


1 4 cay UT 444 _ qT, . UT. . 


~ a cos 6, 5B Ad 


Ty 4M, 2415089.) - Tos 4 42 OP 5-1, 
Ar 


—* ae 7* =k 
al (it) j erie] y inj “Yi-1g | 
3 2AA 


—* mere —* 
Tay um, .,,cos0, ., = Tis uM, ,_,cos8._. 
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The program as written gives the average temperatures, it and ut values 
at velocity points for the second order differencing. For fourth order 
differencing the grid distance is doubled and the values T and un are 
determined at the adjacent mass points. All the averaged quantities 
are determined similarly to those in Appendix A. 

The flux term of the continuity equation is: 

ona OMLik F UM 5k 7 

3 Ar 


a cos 6, 
J 


VT 5 4.608955 _ V4 ink" 5-14 
Ar 
~r* a 
2 thieajk 7 SMinajk | 
3 2Ar 
VT 5 «41, C088. 4 - V4 j-2°08°5-2,, 
2A0 


The averaged values are also determined as indicated in Appendix A. 
The unstarred quantities are at a velocity point and the starred 
quantities are at a mass point. For this eS fourth order 
differencing was blended with the second order differencing at 75° 
latitude. This was done to avoid having to go over the poles to pick 


up the j+l, j+2 and j+3 quantities at the higher latitudes. 
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